#include "MultiVariateInterpolation.h"
////////////////////////////////////////////////////////////////////////////////
MultiVariateInterpolation::MultiVariateInterpolation( )
{
   //this->bounds = boundsSet;
}
////////////////////////////////////////////////////////////////////////////////
MultiVariateInterpolation::~MultiVariateInterpolation( )
{
}
////////////////////////////////////////////////////////////////////////////////
void MultiVariateInterpolation::SetInputPoints( 
      
      map< int, vector< double > > inp, multi_array<double,3> out )
{
   this->inputMap = inp;
   //set shape of outArray
   //array outputArray(extents[bounds.at(0)][bounds.at(1)][bounds.at(2)]);
   array3 outputArray = out;
   for ( unsigned int i=0;i<sizeof( outputArray.shape( ) )-1;i++ )
   {
      if ( DBG_MV_INTERP )
         cout<<" DBG : Size of s "<<i<<":"<<" "<<outputArray.shape( )[i]<<endl;
         bounds.push_back( outputArray.shape( )[i] );
   }
   if ( DBG_MV_INTERP > 2 )
   {
      for ( unsigned int i=0;i<sizeof( outputArray.shape( ) )-1;i++ )
      {
         cout<<" DBG :   BOUNDS : "<<bounds.at(i)<<endl;
      }
   }
   if ( DBG_MV_INTERP )
   {
      for ( unsigned int i=0;i<outputArray.shape()[0];i++ )
      {
         for ( unsigned int j=0;j<outputArray.shape()[1];j++ )
         {
            for ( unsigned int p=0;p<outputArray.shape()[2];p++ )
            {
               cout<<i<<"\t"<<j<<"\t"<<p<<"\t"<<outputArray[i][j][p]<<endl;
            }
         }
      }
   }
   //assume 4 as the dimension for now until I learn more about boost
   //the vector of vectors represents each set of input parameters
   //for the n-dimensional splines
   if ( DBG_MV_INTERP )
   {
      cout<<" DBG : Number of parameters is :"<<inputMap.size()<<endl;
      cout<<" DBG : Number of outputs is :"<<out.num_elements()<<endl;
      cout<<" DBG : Size of multi array is :"<<out.size()<<endl;
   }
   if ( inputMap.size() != out.num_dimensions() )
   {
      cout<<" DBG : Please check input/output data mapping !!!"<<endl;
      exit( -1 );
   }
}
////////////////////////////////////////////////////////////////////////////////
void MultiVariateInterpolation::SetNumberOfParameters( int params )
{
   this->numOfParameters = params;
}
////////////////////////////////////////////////////////////////////////////////

double MultiVariateInterpolation::EvaluateCubicSplineAt( 
   vector<double> evalPoints, multi_array<double,3> out )
   
{

   if ( DBG_MV_INTERP > 3 )
   {
      cout<<" DBG : EVALUATING !!!!!!!!!!!!!"<<endl;
      for ( int i=0;i<bounds.at(0);i++ )
      {
         for ( int j=0;j<bounds.at(1);j++ )
         {
            for ( int p=0;p<bounds.at(2);p++ )
            {
               cout<<" DBG : Here :"<<
                     i<<"\t"<<j<<"\t"<<p<<"\t"<<out[i][j][p]<<endl;
            }
         }
      }
   }

   //set the evaluation points
   evaluationPoints.clear( );
   for ( unsigned int i=0;i<evalPoints.size();i++ )
   {
      evaluationPoints.push_back( evalPoints.at(i) );
   }
   if ( (int) evaluationPoints.size() > numOfParameters )
   {
      
      cout<<" DBG : Evaluation vector has more components than dimensionality"<<
            " of data :  Please check !!!"<<endl;
      exit( -1 );
   }
   else if ( (int) evaluationPoints.size() < numOfParameters )
   {
      cout<<" DBG : Evaluation vector is smaller than dimensionality of data, "
            <<"remaining dimensions will be set to 0"<<endl;
      
      evaluationPoints.resize( numOfParameters );
      for ( int i=evaluationPoints.size();i<numOfParameters;i++ )
      {
         evaluationPoints.at(i) = 0;
      }
   }
   //loop over the number of parameters, which is the same as the number of
   //directions
   vector<double> xdata, ydata;
   vector<double> resultData;
   double tmpRes;
   double fVal;
   fVal = 0;
   //get the xdata and corresponding ydata values to fit 1D splines
   xdata = inputMap.find(0)->second;
   if ( DBG_MV_INTERP > 3 )
   {
      for ( unsigned int i=0;i<xdata.size();i++ )
         
         cout<<" DBG : XDATA : "<<xdata.at(i)<<endl;
   }
   ydata.resize( xdata.size() );
   for ( int k=0;k<bounds.at(2);k++ )
   {
      for ( int j=0;j<bounds.at(1);j++ )
      {
         for ( int i=0;i<bounds.at(0);i++ )
         {
            ydata.at( i ) = out[i][j][k];
            if ( DBG_MV_INTERP > 3 )
               cout<<" DBG :  YDATA : "<<i<<"\t"<<j<<"\t"<<k<<"\t"<<
                     out[i][j][k]<<endl;
         }
         spline1D = new Interpolation( );
         spline1D->SetInputPoints1D( xdata, ydata );
         tmpRes = spline1D->EvaluateCubicSplineAt( evaluationPoints.at(0)
               , xdata );
         
         if ( DBG_MV_INTERP > 3 )
            cout<<" DBG :   TMPRES : "<<tmpRes<<endl;
         resultData.push_back( tmpRes );
         //clean up
         delete spline1D;
         
         //cin >> xat;
      }
   }
   if ( DBG_MV_INTERP > 3 )
      cout<<" DBG :   Size of results : "<<resultData.size( )<<endl;
   if ( DBG_MV_INTERP > 3 )
   {
      for ( unsigned int c=0;c<resultData.size( );c++ )
         cout<<" DBG :   RESULTS DATA : "<<c<<" "<<resultData.at(c)<<endl;
   }
   //splines in the y direction
   xdata.clear( );
   
   xdata = inputMap.find(1)->second;
   ydata.resize( xdata.size() );
   if ( DBG_MV_INTERP > 3 )
      cout<<" DBG :   Size of xdata and ydata (y) : "<<xdata.size( )<<
            " "<<ydata.size()<<endl;
   int counter = 0;
   vector<double> resultsY;
   
   resultsY.clear( );
   for ( unsigned int i=0;i<resultData.size( )+1;i++ )
   {
      if ( counter == (int)xdata.size() )
      {
         counter = 0;
         if ( DBG_MV_INTERP > 3 )
            cout<<" DBG : Fitting spline : "<<i-xdata.size()<<endl;
         spline1D = new Interpolation( );
         spline1D->SetInputPoints1D( xdata, ydata );
         tmpRes = spline1D->
               EvaluateCubicSplineAt( evaluationPoints.at(1), xdata );
         if ( DBG_MV_INTERP > 3 )
            cout<<" DBG :   TMPRES : "<<tmpRes<<endl;
         resultsY.push_back( tmpRes );
         ydata.clear( );
         ydata.resize( xdata.size() );
         //clean up
         delete spline1D;
         if ( DBG_MV_INTERP > 3 )
            cout<<" DBG : Size of resultsY : "<<resultsY.size( )<<endl;
      }
      else
      {
         ydata.at(counter) = resultData.at(i);
      }
      counter++;
   }
   if ( DBG_MV_INTERP > 3 )
      cout<<" DBG : Size of resultsY : Finished for loop : "<<resultsY.size( )
            <<endl;
   
   if ( DBG_MV_INTERP > 2 )
   {
      
      for ( unsigned int c=0;c<resultsY.size( );c++ )
         cout<<" DBG :   RESULTS Y : "<<resultsY.at(c)<<" C : "<<c<<endl;
   }
   xdata.clear( );
   ydata.clear( );
   xdata = inputMap.find(2)->second;
   ydata.resize( xdata.size() );
   if ( DBG_MV_INTERP > 3 )
      
   {
      cout<<" DBG :  Setting data for next spline now :"<<endl;
      cout<<" DBG :  X DATA SIZE : "<<xdata.size( )<<endl;
      cout<<" DBG :  Y DATA SIZE : "<<ydata.size( )<<endl;
   }
   for ( unsigned int i=0;i<xdata.size();i++ )
   {
      ydata.at(i) = resultsY.at(i);
   }
   spline1D = new Interpolation( );
   spline1D->SetInputPoints1D( xdata, ydata );
   tmpRes = spline1D->EvaluateCubicSplineAt( evaluationPoints.at(2), 
         xdata );
   if ( DBG_MV_INTERP > 2 )
      cout<<" DBG : RESULT : "<<tmpRes<<endl;
   //clean up
   delete spline1D;

   //clean up
   xdata.clear( );
   ydata.clear( );
   
   return tmpRes;
}

////////////////////////////////////////////////////////////////////////////////
void MultiVariateInterpolation::SetInputPoints( 
      
      map< int, vector< double > > inp, multi_array<double,4> out )
{
   this->inputMap = inp;
   //set shape of outArray
   array4 outputArray = out;
   for ( unsigned int i=0;i<sizeof( outputArray.shape( ) );i++ )
   {
      if ( DBG_MV_INTERP )
         cout<<" DBG : Size of s "<<i<<":"<<" "<<outputArray.shape( )[i]<<endl;
         bounds.push_back( outputArray.shape( )[i] );
   }
   if ( DBG_MV_INTERP > 2 )
   {
      cout<<" DBG : SIZE OF BOUNDS : "<<bounds.size( )<<endl;
      for ( unsigned int i=0;i<sizeof( outputArray.shape( ) );i++ )
      {
         cout<<" DBG :   BOUNDS : "<<i<<" "<<bounds.at(i)<<endl;
      }
   }
   if ( DBG_MV_INTERP )
   {
      for ( unsigned int i=0;i<outputArray.shape()[0];i++ )
      {
         for ( unsigned int j=0;j<outputArray.shape()[1];j++ )
         {
            for ( unsigned int p=0;p<outputArray.shape()[2];p++ )
            {
               for ( unsigned int q=0;q<outputArray.shape()[3];q++ )
               {
                  cout<<i<<"\t"<<j<<"\t"<<p<<"\t"<<q<<"\t"<<
                        outputArray[i][j][p][q]<<endl;
               }
            }
         }
      }
   }
   //assume 4 as the dimension for now until I learn more about boost
   //the vector of vectors represents each set of input parameters
   //for the n-dimensional splines
   if ( DBG_MV_INTERP )
   {
      cout<<" DBG : Number of parameters is :"<<inputMap.size()<<endl;
      cout<<" DBG : Number of outputs is :"<<out.num_elements()<<endl;
      cout<<" DBG : Size of multi array is :"<<out.size()<<endl;
   }
   if ( inputMap.size() != out.num_dimensions() )
   {
      cout<<" DBG : Please check input/output data mapping !!!"<<endl;
      exit( -1 );
   }
}
////////////////////////////////////////////////////////////////////////////////
double MultiVariateInterpolation::EvaluateCubicSplineAt
   
   ( vector<double> evalPoints, multi_array<double,4> out )
{
   if ( DBG_MV_INTERP > 3 )
   {
      cout<<" DBG : EVALUATING !!!!!!!!!!!!!"<<endl;
      for ( int i=0;i<bounds.at(0);i++ )
      {
         for ( int j=0;j<bounds.at(1);j++ )
         {
            for ( int p=0;p<bounds.at(2);p++ )
            {
               for ( int q=0;q<bounds.at(3);q++ )
                  cout<<" DBG : Here :"<<
                     i<<"\t"<<j<<"\t"<<p<<"\t"<<q<<"\t"<<out[i][j][p][q]<<endl;
            }
         }
      }
   }

   
   //set the evaluation points
   evaluationPoints.clear( );
   
   for ( unsigned int i=0;i<evalPoints.size();i++ )
   {
      evaluationPoints.push_back( evalPoints.at(i) );
   }
   if ( DBG_MV_INTERP > 2 )
      cout<<" DBG : Number of evaluation points : "<<
            evaluationPoints.size( )<<endl;
   
   if ( (int)evaluationPoints.size() > numOfParameters )
   {
      cout<<" DBG : Evaluation vector has more components than dimensionality"<<
            " of data  Please check !!!"<<endl;
      exit( -1 );
   }
   else if ( (int)evaluationPoints.size() < numOfParameters )
   {
      cout<<" DBG : Evaluation vector is smaller than dimensionality of data, "
            <<"remaining dimensions will be set to 0"<<endl;
      evaluationPoints.resize( numOfParameters );
      for ( int i=evaluationPoints.size();i<numOfParameters;i++ )
      {
         evaluationPoints.at(i) = 0;
      }
   }
   //let the directions be x1, x2, x3 and x4
   vector<double> xdata, ydata;
   vector<double> result1Data, result2Data;
   double tmpRes;
   xdata = inputMap.find(0)->second;
   if ( DBG_MV_INTERP > 3 )
   {
      for ( unsigned int i=0;i<xdata.size();i++ )
         
         cout<<" DBG : XDATA : "<<xdata.at(i)<<endl;
   }
   ydata.resize( xdata.size() );
   for ( int l=0;l<bounds.at(3);l++ )
   {
      for ( int k=0;k<bounds.at(2);k++ )
      {
         for ( int j=0;j<bounds.at(1);j++ )
         {
            for ( int i=0;i<bounds.at(0);i++ )
            {
               ydata.at( i ) = out[i][j][k][l];
               if ( DBG_MV_INTERP > 3 )
                  cout<<" DBG :  YDATA : "<<i<<"\t"<<j<<"\t"<<k<<"\t"<<l<<"\t"<<
                     out[i][j][k][l]<<endl;
            }
            spline1D = new Interpolation( );
            spline1D->SetInputPoints1D( xdata, ydata );
            tmpRes = spline1D->EvaluateCubicSplineAt( evaluationPoints.at(0)
               , xdata );
            
            if ( DBG_MV_INTERP > 3 )
               cout<<" DBG :   TMPRES : "<<tmpRes<<endl;
            result1Data.push_back( tmpRes );
            //clean up
            
            delete spline1D;
            //cin >> xat;
         }
      }
   }
   if ( DBG_MV_INTERP > 2 )
      cout<<" DBG :   Size of results : "<<result1Data.size( )<<endl;
   if ( DBG_MV_INTERP > 3 )
   {
      
      for ( unsigned int c=0;c<result1Data.size( );c++ )
         cout<<" DBG :   RESULTS DATA : "<<c<<" "<<result1Data.at(c)<<endl;
   }
   //fit splines in the x2 direction
   xdata.clear( );
   xdata = inputMap.find(1)->second;
   ydata.clear( );
   
   ydata.resize( xdata.size() );
   //this is a 3D variables problem now, create a input/output map and 
   //call the 3D function
   array3 out3D(extents[bounds.at(1)][bounds.at(2)][bounds.at(3)]);
   int counter = 0;
   for ( int k=0;k<bounds.at(3);k++ )
   {
      for ( int j=0;j<bounds.at(2);j++ )
      {
         for ( int i=0;i<bounds.at(1);i++ )
         {
            out3D[i][j][k] = result1Data.at(counter);
            counter++;
         }
      }
   }
   if ( DBG_MV_INTERP > 2 )
      
      cout<<" DBG : Counter : "<<counter<<endl;
   //create 3D map
   map<int, vector<double> > map3D;
   vector<double> x2, x3, x4;
   vector<double> eval3D;
   eval3D.push_back( evaluationPoints.at(1) );
   eval3D.push_back( evaluationPoints.at(2) );
   eval3D.push_back( evaluationPoints.at(3) );
   
   if ( DBG_MV_INTERP )
   {
      cout<<" DBG : Created map "<<endl;
      cout<<" DBG : Size of evaluation vector : "<<eval3D.size( )<<endl;
   }
   x2 = inputMap.find(1)->second;
   x3 = inputMap.find(2)->second;
   x4 = inputMap.find(3)->second;
   map3D.insert( make_pair( 0, x2 ) );
   map3D.insert( make_pair( 1, x3 ) );
   map3D.insert( make_pair( 2, x4 ) );
   
   
   MultiVariateInterpolation* mv3 = new MultiVariateInterpolation( );
   mv3->SetNumberOfParameters( 3 );
   mv3->SetInputPoints( map3D, out3D );
   
   tmpRes = mv3->EvaluateCubicSplineAt( eval3D, out3D );
   
   return tmpRes;
}

////////////////////////////////////////////////////////////////////////////////
double MultiVariateInterpolation::EvaluateLinearAt( vector<double> evalPoints, 
         multi_array<double,3> out )
{
   //fit a line in x1 direction
   
      if ( DBG_MV_INTERP > 3 )
   {
      cout<<" DBG : EVALUATING !!!!!!!!!!!!!"<<endl;
      for ( int i=0;i<bounds.at(0);i++ )
      {
         for ( int j=0;j<bounds.at(1);j++ )
         {
            for ( int p=0;p<bounds.at(2);p++ )
            {
               cout<<" DBG : Here :"<<
                     i<<"\t"<<j<<"\t"<<p<<"\t"<<out[i][j][p]<<endl;
            }
         }
      }
   }

   //set the evaluation points
   evaluationPoints.clear( );
   for ( unsigned int i=0;i<evalPoints.size();i++ )
   {
      evaluationPoints.push_back( evalPoints.at(i) );
   }
   if ( (int)evaluationPoints.size() > numOfParameters )
   {
      
      cout<<" DBG : Evaluation vector has more components than dimensionality"<<
            " of data :  Please check !!!"<<endl;
      exit( -1 );
   }
   else if ( (int)evaluationPoints.size() < numOfParameters )
   {
      cout<<" DBG : Evaluation vector is smaller than dimensionality of data, "
            <<"remaining dimensions will be set to 0"<<endl;
      
      evaluationPoints.resize( numOfParameters );
      for ( int i=evaluationPoints.size();i<numOfParameters;i++ )
      {
         evaluationPoints.at(i) = 0;
      }
   }
   //loop over the number of parameters, which is the same as the number of
   //directions
   vector<double> xdata, ydata;
   vector<double> resultData;
   double tmpRes;
   double fVal;
   fVal = 0;
   //get the xdata and corresponding ydata values to fit 1D splines
   xdata = inputMap.find(0)->second;
   if ( DBG_MV_INTERP > 3 )
   {
      for ( unsigned int i=0;i<xdata.size();i++ )
         
         cout<<" DBG : XDATA : "<<xdata.at(i)<<endl;
   }
   ydata.resize( xdata.size() );
   for ( int k=0;k<bounds.at(2);k++ )
   {
      for ( int j=0;j<bounds.at(1);j++ )
      {
         for ( int i=0;i<bounds.at(0);i++ )
         {
            ydata.at( i ) = out[i][j][k];
            if ( DBG_MV_INTERP > 3 )
               cout<<" DBG :  YDATA : "<<i<<"\t"<<j<<"\t"<<k<<"\t"<<
                     out[i][j][k]<<endl;
         }
         line1D = new Interpolation( );
         line1D->SetInputPoints1D( xdata, ydata );
         tmpRes = line1D->EvaluateLinearAt( evaluationPoints.at(0), xdata );
         
         
         if ( DBG_MV_INTERP > 3 )
            cout<<" DBG :   TMPRES : "<<tmpRes<<endl;
         resultData.push_back( tmpRes );
         //clean up
         delete line1D;
         
         //cin >> xat;
      }
   }
   if ( DBG_MV_INTERP > 3 )
      cout<<" DBG :   Size of results : "<<resultData.size( )<<endl;
   if ( DBG_MV_INTERP > 3 )
   {
      for ( unsigned int c=0;c<resultData.size( );c++ )
         cout<<" DBG :   RESULTS DATA : "<<c<<" "<<resultData.at(c)<<endl;
   }
   //splines in the y direction
   xdata.clear( );
   
   xdata = inputMap.find(1)->second;
   ydata.resize( xdata.size() );
   if ( DBG_MV_INTERP > 3 )
      cout<<" DBG :   Size of xdata and ydata (y) : "<<xdata.size( )<<
            " "<<ydata.size()<<endl;
   int counter = 0;
   vector<double> resultsY;
   
   resultsY.clear( );
   for ( unsigned int i=0;i<resultData.size( )+1;i++ )
   {
      if ( counter == (int)xdata.size() )
      {
         counter = 0;
         if ( DBG_MV_INTERP > 3 )
            cout<<" DBG : Fitting spline : "<<i-xdata.size()<<endl;
         line1D = new Interpolation( );
         line1D->SetInputPoints1D( xdata, ydata );
         tmpRes = line1D->EvaluateLinearAt( evaluationPoints.at(1), xdata );
         if ( DBG_MV_INTERP > 3 )
            cout<<" DBG :   TMPRES : "<<tmpRes<<endl;
         resultsY.push_back( tmpRes );
         ydata.clear( );
         ydata.resize( xdata.size() );
         //clean up
         delete line1D;
         if ( DBG_MV_INTERP > 3 )
            cout<<" DBG : Size of resultsY : "<<resultsY.size( )<<endl;
      }
      else
      {
         ydata.at(counter) = resultData.at(i);
      }
      counter++;
   }
   if ( DBG_MV_INTERP > 3 )
      cout<<" DBG : Size of resultsY : Finished for loop : "<<resultsY.size( )
            <<endl;
   
   if ( DBG_MV_INTERP > 2 )
   {
      for ( unsigned int c=0;c<resultsY.size( );c++ )
         cout<<" DBG :   RESULTS Y : "<<resultsY.at(c)<<" C : "<<c<<endl;
   }
   xdata.clear( );
   ydata.clear( );
   xdata = inputMap.find(2)->second;
   ydata.resize( xdata.size() );
   if ( DBG_MV_INTERP > 3 )
      
   {
      cout<<" DBG :  Setting data for next spline now :"<<endl;
      cout<<" DBG :  X DATA SIZE : "<<xdata.size( )<<endl;
      cout<<" DBG :  Y DATA SIZE : "<<ydata.size( )<<endl;
   }
   for ( unsigned int i=0;i<xdata.size();i++ )
   {
      ydata.at(i) = resultsY.at(i);
   }
   line1D = new Interpolation( );
   line1D->SetInputPoints1D( xdata, ydata );
   tmpRes = line1D->EvaluateLinearAt( evaluationPoints.at(2), 
         xdata );
   if ( DBG_MV_INTERP > 2 )
      cout<<" DBG : RESULT : "<<tmpRes<<endl;
   //clean up
   delete line1D;

   //clean up
   xdata.clear( );
   ydata.clear( );
   
   return tmpRes;
   
}
////////////////////////////////////////////////////////////////////////////////
double MultiVariateInterpolation::EvaluateLinearAt( vector<double> evalPoints, 
   
   multi_array<double,4> out )
{
   
   
      if ( DBG_MV_INTERP > 3 )
   {
      cout<<" DBG : EVALUATING !!!!!!!!!!!!!"<<endl;
      for ( int i=0;i<bounds.at(0);i++ )
      {
         for ( int j=0;j<bounds.at(1);j++ )
         {
            for ( int p=0;p<bounds.at(2);p++ )
            {
               for ( int q=0;q<bounds.at(3);q++ )
                  cout<<" DBG : Here :"<<
                     i<<"\t"<<j<<"\t"<<p<<"\t"<<q<<"\t"<<out[i][j][p][q]<<endl;
            }
         }
      }
   }

   
   //set the evaluation points
   evaluationPoints.clear( );
   
   for ( unsigned int i=0;i<evalPoints.size();i++ )
   {
      evaluationPoints.push_back( evalPoints.at(i) );
   }
   if ( DBG_MV_INTERP > 2 )
      cout<<" DBG : Number of evaluation points : "<<
            evaluationPoints.size( )<<endl;
   
   if ( (int)evaluationPoints.size() > numOfParameters )
   {
      cout<<" DBG : Evaluation vector has more components than dimensionality"<<
            " of data  Please check !!!"<<endl;
      exit( -1 );
   }
   else if ( (int)evaluationPoints.size() < numOfParameters )
   {
      cout<<" DBG : Evaluation vector is smaller than dimensionality of data, "
            <<"remaining dimensions will be set to 0"<<endl;
      evaluationPoints.resize( numOfParameters );
      for ( int i=evaluationPoints.size();i<numOfParameters;i++ )
      {
         evaluationPoints.at(i) = 0;
      }
   }
   //let the directions be x1, x2, x3 and x4
   vector<double> xdata, ydata;
   vector<double> result1Data, result2Data;
   double tmpRes;
   xdata = inputMap.find(0)->second;
   if ( DBG_MV_INTERP > 3 )
   {
      for ( unsigned int i=0;i<xdata.size();i++ )
         
         cout<<" DBG : XDATA : "<<xdata.at(i)<<endl;
   }
   ydata.resize( xdata.size() );
   for ( int l=0;l<bounds.at(3);l++ )
   {
      for ( int k=0;k<bounds.at(2);k++ )
      {
         for ( int j=0;j<bounds.at(1);j++ )
         {
            for ( int i=0;i<bounds.at(0);i++ )
            {
               ydata.at( i ) = out[i][j][k][l];
               if ( DBG_MV_INTERP > 3 )
                  cout<<" DBG :  YDATA : "<<i<<"\t"<<j<<"\t"<<k<<"\t"<<l<<"\t"<<
                     out[i][j][k][l]<<endl;
            }
            line1D = new Interpolation( );
            line1D->SetInputPoints1D( xdata, ydata );
            tmpRes = line1D->EvaluateLinearAt( evaluationPoints.at(0), xdata );
            
            if ( DBG_MV_INTERP > 3 )
               cout<<" DBG :   TMPRES : "<<tmpRes<<endl;
            result1Data.push_back( tmpRes );
            //clean up
            
            delete line1D;
            //cin >> xat;
         }
      }
   }
   if ( DBG_MV_INTERP > 2 )
      cout<<" DBG :   Size of results : "<<result1Data.size( )<<endl;
   if ( DBG_MV_INTERP > 3 )
   {
      
      for ( unsigned int c=0;c<result1Data.size( );c++ )
         cout<<" DBG :   RESULTS DATA : "<<c<<" "<<result1Data.at(c)<<endl;
   }
   //fit splines in the x2 direction
   xdata.clear( );
   xdata = inputMap.find(1)->second;
   ydata.clear( );
   
   ydata.resize( xdata.size() );
   //this is a 3D variables problem now, create a input/output map and 
   //call the 3D function
   array3 out3D(extents[bounds.at(1)][bounds.at(2)][bounds.at(3)]);
   int counter = 0;
   for ( int k=0;k<bounds.at(3);k++ )
   {
      for ( int j=0;j<bounds.at(2);j++ )
      {
         for ( int i=0;i<bounds.at(1);i++ )
         {
            out3D[i][j][k] = result1Data.at(counter);
            counter++;
         }
      }
   }
   if ( DBG_MV_INTERP > 2 )
      
      cout<<" DBG : Counter : "<<counter<<endl;
   //create 3D map
   map<int, vector<double> > map3D;
   vector<double> x2, x3, x4;
   vector<double> eval3D;
   eval3D.push_back( evaluationPoints.at(1) );
   eval3D.push_back( evaluationPoints.at(2) );
   eval3D.push_back( evaluationPoints.at(3) );
   
   if ( DBG_MV_INTERP )
   {
      cout<<" DBG : Created map "<<endl;
      cout<<" DBG : Size of evaluation vector : "<<eval3D.size( )<<endl;
   }
   x2 = inputMap.find(1)->second;
   x3 = inputMap.find(2)->second;
   x4 = inputMap.find(3)->second;
   map3D.insert( make_pair( 0, x2 ) );
   map3D.insert( make_pair( 1, x3 ) );
   map3D.insert( make_pair( 2, x4 ) );
   
   
   MultiVariateInterpolation* mv3 = new MultiVariateInterpolation( );
   mv3->SetNumberOfParameters( 3 );
   mv3->SetInputPoints( map3D, out3D );
   tmpRes = mv3->EvaluateLinearAt( eval3D, out3D );
   
   return tmpRes;
}
